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We present a method for reconstructing the photon num- 
ber distribution from the homodyne statistics based on maxi- 
mization of the likelihood function derived from the exact sta- 
tistical description of a homodyne experiment. This method 
incorporates in a natural way the physical constraints on 
the reconstructed quantities, and the compensation for the 
nonunit detection efficiency 

PACS Number(s): 42.50.Ar, 42.50.Dv 



An interesting application of the homodyne detection 
is the measurement of phase-insensitive properties of op- 
tical fields, dependent on the photon number distribution 
0. The homodyne technique goes beyond the current 
limitations of direct photodetection. First, ultrafast sam- 
pling time can be achieved by using the local oscillator 
field in the form of a short pulse. Secondly, information 
on the photon distribution is carried by two relatively 
intense fields, which can be detected with substantially 
higher efficiency than the signal field itself. This fea- 
ture has enabled an experimental demonstration of even- 
odd oscillations in the photon number distribution of a 
squeezed vacuum state Q . 

In the homodyne scheme, all phase-independent prop- 
erties of the measured field are contained in the phase- 
averaged statistics of difference counts [p). The proba- 
bility distribution of difference counts is a linear combi- 
nation of diagonal elements of the density matrix in the 
Fock basis. This relation can be analytically inverted 
which yields an expression for the photon number 
distribution as integrals of the homodyne statistics with 
the so-called pattern functions ||. In a real experiment, 
however, the homodyne statistics is known only with 
some statistical uncertainty, as it is determined from a 
finite number of experimental runs. Application of pat- 
tern functions to noisy experimental data can generate 
unphysical results, such as negativities in the photon 
number distribution. These artifacts become particularly 
strong, when compensation for the detector imperfection 
is used in the numerical processing of the experimental 
data §]. 

In this communication we apply the statistical method- 
ology of the maximum likelihood estimation |7j|| to re- 
construct the photon number distribution from the ho- 
modyne statistics. This approach incorporates in a natu- 
ral way the finite character of the experimental data sam- 
ple, as well as physical constraints on the reconstructed 
quantities. Furthermore, compensation for the detector 



imperfection is consistently built into the reconstruction 
scheme, without generating unphysical artifacts. Com- 
pared to the recent application of the least-square inver- 
sion method to quantum state reconstruction our 
algorithm is based directly on the exact statistical anal- 
ysis of a homodyne experiment. This automatically as- 
signs proper statistical weights to experimental frequen- 
cies of homodyne events, and no simplifying assumptions 
about the noise present in the finite sample of data are 
necessary. 

We will start with a statistical description of data col- 
lected in a homodyne setup. The phase-averaged proba- 
bility distribution p(x) of recording a specific outcome x 
is a linear combination of the photon number distribution 
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with coefficients given by the formula |10|: 
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where r\ is the detection efficiency and H n (x) denote Her- 
mite polynomials. The continuous variable x is divided 
into bins of the width Ax, which we will label with v. 
When bins are sufficiently small, we may approximate 
the probability p v of registering the outcome in a v\h bin 
by: 



p v = p(x v )Ax, 



(3) 



where x v is the central point of the vth bin. 

Repeating the measurement N times yields a frequency 
histogram {k u } specifying in how many runs the outcome 
has been found in a vth. bin. The probability of obtaining 
a specific histogram {k„} is given by the multinomial 
distribution: 



n{k»}\{Pn})=N\l[ 
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where we have explicitly written its dependence on the 
photon number distribution {p n } entering the right hand 
side via probabilities p v . 

When processing the experimental data, we face an in- 
verse problem: given a certain histogram {k v } we want 
to reconstruct the photon number distribution {p n }- The 
answer given to this problem by the maximum likelihood 
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method is that the best estimate for {p n } maximizes the 
function defined in Eq. (^), with k v 's treated as fixed pa- 
rameters obtained from an experiment. The search for 
the maximum of the function V, called the likelihood 
function, is a priori restricted to the manifold of {p n } 
that describe a possible physical situation. This guaran- 
tees that the reconstructed probability distribution will 
be free from unphysical artifacts. 

This maximization problem has to be solved by nu- 
merical means. For this purpose we will introduce a cut- 
off parameter K for the photon number. The positivity 
constraints for p n can be satisfied by a substitution of 
variables: p n — y^- Instead of computing the likelihood 
function, it is more convenient to consider its logarithm: 



£(2/0)2/1, • • • , Vk) = E K log ^2 A vnyl 
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where we have omitted terms independent of y n , and de- 
noted A un — A n (x u )l±.x. The condition that the sum of 
probabilities p n is equal to one can be taken into account 
using a Lagrange multiplier A/\ Thus we obtain a set of 
K + 1 equations: 



can be applied to incomplete data consisting only of se- 
lected histogram bins, provided that they contain enough 
information to resolve contributions from different Fock 
states. This feature may be useful when, for example, 
statistics in some bins is corrupted by external noise. 

In conclusion, we have presented a method for recon- 
structing the photon number distribution from homo- 
dyne data via maximization of the likelihood function. 
It allows one to reduce the statistical error by includ- 
ing in the reconstruction scheme a priori constraints on 
the quantities to be determined. This method has solid 
methodological background and has been derived from 
the exact statistical description of a homodyne experi- 
ment. 

The author is indebted to Professor K. Wodkiewicz 
for numerous discussions and valuable comments on the 
manuscript. This research was supported by the Polish 
KBN grant No. 2P03B 002 14 and by the Foundation for 
Polish Science. 
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Multiplying these equations by y m and adding them to- 
gether yields: 
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that is, the Lagrange multiplier is equal to the total num- 
ber of experimental runs. 

The maximization problem reformulated in this way 
can be treated with standard numerical procedures. We 
have performed Monte Carlo simulations of the ho- 
modyne experiment and applied the downhill simplex 
method to reconstruct the photon number distri- 
bution via maximization of the likelihood function. In 
Fig. |l| we depict estimation of the photon number distri- 
bution for a coherent state and a squeezed vacuum state, 
both states with the average photon number equal to 
one. We have assumed imperfect detection characterized 
by the efficiency r\ = 85%, which has been taken into 
account in the reconstruction process by setting appro- 
priately coefficients A n (x v ) defined in Eq. (||). The sim- 
ulated homodyne data are depicted in the top graphs, 
and the result of the reconstruction is compared with 
theoretical photon number distributions in the bottom 
graphs. The even-odd oscillations in the photon num- 
ber distribution of the squeezed vacuum state are fully 
recovered despite losses in the detection process. Let us 
note that the maximum likelihood estimation algorithm 
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FIG. 1. Monte Carlo simulations for a coherent state (left) 
and a squeezed vacuum state (right). The top graphs depict 
phase-averaged homodyne statistics obtained from N = 10 5 
runs, superposed on analytical distributions shown with solid 
lines. The interval —5 < x < 5 is divided into 100 bins. 
The bottom graphs compare the reconstructed photon num- 
ber distribution (circles) with analytical values (bars). The 
cut-off parameter is K — 19. 
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